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Computation of propagation effects in the neutral atmosphere, namely path delay, extinction, and 
bending angle is a trivial task provided the 4D state of the atmosphere is known. Unfortunately, 
the mixing ratio of water vapor is highly variable and it cannot be deduced from surface measure¬ 
ments. That fact led to a paradigm that considers path delay and extinction in the atmosphere as 
a priori unknown quantities that have to be evaluated from the radio astronomy data themselves. 
Development of our ability to model the atmosphere and to digest humongous outputs of these 
models that took place over the course of the 21st century changed the game. Using the publicly 
available output of operational numerical weather model GEOS run by NASA, we are in a po¬ 
sition to compute path delay through the neutral atmosphere for any station and for any epoch 
from 1979 through now with accuracy of 45 ps * cosec elevation. We are in a position to compute 
extinction with accuracy better than 10 pro cents. We are in a position to do it routinely, in a 
similar way how we update apparent star positions for precession and nutation. Moreover, we are 
in a position to do it now. As a demonstration of current capabilities, I have computed time series 
of path delays for aall radiotelecopes that I was aware of (220 sites) since 1979 with a step 3-6 
hours. Results of the validation tests are presented. A new paradigm of data analysis assumes 
that we know the atmosphere propagation effects a priori with the accuracy higher that one could 
deduce them from radio astronomy observations. 
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1. Introduction 

Radio waves travel billions years from typical sources observed with VLBI. At the very end of 
their journey, at the last millisecond, they propagate through the Earth’s atmosphere. While multi¬ 
frequency observations allow us to evaluate the contribution of the ionized component of the Earth’s 
atmosphere, modeling propagation effects in the neutral atmosphere poses a challenge. There are 
four effects: 1) refraction, i.e. trajectory bending, 2) delay in the atmosphere, 3) attenuation, and 
4) atmosphere emission. If we know the state of the atmosphere, we can compute these quantities. 
The crux of the problem is that the state of the 4D atmosphere cannot be deduced from surface 
measurements. The reason of that is water is underwent phase transition at the atmospheric layers 
that are above the surface, and as a result, the mixing ratio of water vapor is highly volatile. 

2. Old paradigm 

The old paradigm assumed we do not know the instantaneous state of the atmosphere. In 
order to circumvent the lack of knowledge, some greatly simplified or just wrong a priori models 
were used for computation of path delay in the past. In particular, Saastamoinen (1972) model 
became very popular. That model, in a way how it was implemented, assumed that water vapor 
present in the atmosphere had thermodynamic properties of dry air. A typical error of path delay 
computed that way was 300-1300 ps, i.e. 5-15% of the total delay. In order to improve modeling, 
the astronomical data themselves were suggested to be used for solving residual parameters of the 
propagation model. It is assumed in the framework of this paradigm that the path delay z(t,e,A ) 
can be 1 decomposed into two azimuth-independent components: 

T (t,e,A) = Zd(t)-m d (e) + t w (t)-m w (e), (2.1) 

where T j(t) is the so-called dry component, is the so-called dry mapping function, i.e. a 

derivative with respect to elevation; t w (t) and m w (e) is wet path delay and its derivative with respect 
to elevation angle; t, e, and A stands for time, elevation angle above the horizon, and azimuth. 

For evaluation of extinction in the atmosphere measurements of antenna brightness tempera¬ 
ture at different elevations were often made (tipping curves). 

3. New paradigm 

Advances in numerical weather models made it possible to evaluate parameters of the atmo¬ 
sphere using various ground, air-born, and space-born measurements that are assimilated into a 
dynamic model of the atmosphere. The output of these models defines the parameters of the state 
of the atmosphere on a 4-dimensional grid. Three parameters are important for reduction of astron¬ 
omy observations: air temperature T, total atmospheric pressure P, and partial pressure of water 
vapor P w . Thus, the state of the atmosphere in the new paradigm is considered known. At the 
moment, there are several centers in the world that produce numerical models of the atmosphere. I 
used models produced by the NASA Global Modeling and Assimilation Office (GMAO) 2 for this 
study. 

1 Unfortunately, it is often forgotten, that such a decomposition is a only a simplified approximation within its range 
of applicability. 

“The output of the models is available at http://gmao.gsfc.nasa.gov 
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Table 

1: GMAO models of the atmosphere 

MERRA: 

Since 1979.01.01 

72 lev x0.5° x 0.67° x 6 h 

Latency: 40 J 

GEOS FPIT 

Since 2000.01.01 

72 lev x0.5° x 0.67° x 3 h 

Latency: 1 2 h 

GEOS FP 

Since 2011.09.01 

72 lev x0.25°x 0.3125° x3 ,! 

Latency: I2 h 


4. Computation of path delay 


Electromagnetic waves in the medium propagate with group velocity v m that is always smaller 
than than the speed of light in vacuum c. Propagating through the medium that has variable re- 
fractivity defined as r = the electromagnetic wave changes its direction. Refractivity is a 
simple function of P, T, and P w that depends on coefficients measured in laboratory. According to 
Aparicio & Laroche (2011), refractivity in radio range can be computed with accuracy ~0.1%. 

There are several ways to solve for ray trajectory in a continuous medium with known refrac¬ 
tivity field. I prefer an approach based on solving the variational problem: the actual trajectory that 
the electromagnetic wave travels is the one that minimizes travel time (Fermat principle). This prin¬ 
ciple was first formulated empirically in 1662 but nowadays can be derived directly from Maxwell 
equations (see, for instance Landau and Lifshitz, 1988). A general solution of this problem of cal¬ 
culus of variations was found by Euler in 1744 and it is reduced to solving a non-linear system 
of differential equations of the 4th order with mixed initial values. The coefficients of equations 
depend on the refractivity and its gradient. 

Considering propagation of electromagnetic wave through a continuous medium requires rep¬ 
resentation of refractivity r(h,X,(p g d,t) as a continuous differentiable function. The mathematical 
model that I used for such a representation is an expansion of refractivity into a 4D tensor product 
of B-spline functions of the mth degree fi'": 

r(h,X,<t sd , <)='"£ ‘ £ “IT 'f,iuBT(h)BJ(X)Sp 9s<l )Bn,). (4.1) 

i—l—m 7 = 1 — m k—l—m 1=1—m 


The property of B-spline, minimal support, reduces the problem of finding the coefficients 
fijki to solving a system of algebraic equations with an m-diagonal matrix, which can be solved 
extremely efficiently at modern computers. 

Representation of the refractivity fields in the form of the tensor product of B-splines of the 
3rd degree (m=3), reduces differential equations that are a solution of the variational problem 
to a system of tri-diagonal algebraic non-linear equations. That system is efficiently solved by 
iterations starting from the 1st approximation that radio waves propagate along the straight line. 
Two iterations are sufficient to reduce errors of a numerical solution below 1 ps. 

After determining the trajectory of the radio wave from the emitter to the receiver, we can 
easily find the time delay in the neutral atmosphere along the trajectory tj(£),£(£) by integrating 
refractivity along the path: 
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where the Cartesian coordinate system £,T7,£ is chosen in such a way that axis t, is along the 
straight line from the receiver to the emitter, and axes r \, £ are orthogonal to that direction. 


Attenuation in the atmosphere can be computed by evaluating the integral that resembles [T2 
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(4.3) 


where a(P,P w ,T,f) is the specific attenuation coefficient that has a strong dependency on fre¬ 
quency /. In particular, attenuation of the signal at the ground is a(0). 

Finally, atmosphere brightness temperature T atm is found by solving the radiative transfer equa¬ 
tion, which is reduced to evaluation of the integral 


poo 

Team = / T(Z ,1 7,0 
Jo 





(4.4) 


Figure 1: Baseline length repeatability for a case when residual zenith path delay is solved for (Left) and 
when no path delay was estimated (Right). The baseline length repeatability extrapolated to the Earth’s 
diameter is 1.12 cm for the first case and 2.32 cm for the second case. 




At a given station, given epoch, integrals T na , a, and T atm are functions of azimuth and eleva¬ 
tion. For logistical reasons it is convenient to compute these quantities on a 2D azimuth-elevation 
regular grid (that does not have to be uniform), expand them into the B-spline basis and store the 
coefficients. This approach disentangles computation of integrals EH3 from radio astronomy 
data analysis. 

5. Results 

I computed azimuth-elevation coefficients for path delay expansions for all 220 VLBI stations 
that participated in radio astronomy observations since 1979 through present using GMAO numer¬ 
ical weather models GEOS-FPIT and MERRA. The input dataset has rather significant size, 36Tb, 
but still manageable. 
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In order to evaluate the accuracy of path delays computation, I ran two solutions: with esti¬ 
mating residual path delay in zenith direction and without estimating. For each baseline I estimated 
the weighted root mean square (wrms) of deviations of baseline lengths with respect to the linear 
model — the so-called baseline length repeatability test. The baseline length repeatability depen¬ 
dence on baseline length L (blue lines in Figure [l]) was fitted to R( L) = \JA 1 + (B ■ L) 1 , where 
A and B are fitted coefficients. When the baseline length between two stations is approaching to 
the Earth’s diameter, its vector is approaching to the local vertical at both stations. Therefore, the 
wrms of station position uncertainties can be evaluated as R(L = diam)/\/2. Comparing the base¬ 
line length repeatabilities extrapolated to the Earth’s diameter from two solutions, I derived the 
additional variance in vertical station positions due to errors in a priroi path delays derived from 
the output of numerical weather models: 1.45 cm (48 ps). 

Another approach to evaluate errors of a priori path delay is to compute the rms of estimates of 
residual path delay in zenith direction. It varies from 0.6 cm at polar station NYALES20 to 2.2 cm 
at tropical station SC-VLBA with an average value of 1.25 cm (37 ps). It is instructive to note that 
the average zenith total path delay is 251.4 cm and the average zenith wet path delay is 12.3 cm. 
This allows me to conclude that total path delay can be computed with accuracy 0.5%, and the 
contribution of water vapor to path delay can be computed with accuracy 10%. 

Figure 2: Atmosphere brightness temperature at elevation 45° at SARDINIA at 86.304 GHz 



Feb Apr Jun 


Aug Oct Dec 

Time in 2013 


Figure [2] shows the estimates of the atmosphere brightness temperature at SARDINIA at 
86.304 GHz. This figure may help us to make a rough prediction of the atmosphere contribu¬ 
tion to station performance. We see that a naive surmise that the performance in winter will be 
always significantly better than in summer is a gross simplification. 

Using estimates of the atmosphere attenuation from the output of numerical weather models is 
straightforward: we just calibrate fringe amplitudes by dividing them by e a factor and thus relate 
them to the top of the atmosphere. 
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There are also certain benefits of being able to compute atmosphere brightness temperature 
using the output of numerical weather model. The difference between system temperature T sys and 
T atm is the sum of receiver temperature T rec and the spill-over term that does not depend on time, 
but depends on elevation and, in a less extent, on azimuth. If the receiver works properly, T rec is 
stable. Stacking T sys and T aIm over the time range when T rec is stable, we can separate the spillover 
term from T rec and develop an empirical model for it. Analysis of T rec time series allows us to 
identify the period of time when it had anomalies (see Figure [|). This allows us to clean the data 
for bad estimates of T sys caused by radio interference and restore missing T sys by extrapolating T rec . 


Figure 3: K-band system temperature in K as a function of time in sec (Left) and receiver temperature 
(Right) after subtraction of atmosphere brightness temperature computed using the output of numerical 
weather model GEOS-FPIT. 


Tsys OV-VLBA 



Tsys OV-VLBA 



It is difficult to estimate errors of atmosphere brightness temperature and attenuation in the 
atmosphere directly. It is reasonable to expect that in the frequency range where the dominating 
contribution to the atmosphere attenuation is water vapor, their errors are close to the uncertainty of 
wet path delay, i.e. around 10%. Estimates of accuracies demonstrate that certain old-fashion steps 
of data calibration, such as inserting geodetic blocks into schedules and measuring tipping curves, 
are nowadays unnecessary, since even better accuracies can be achieved by utilizing the output of 
numerical weather models. 

References 


Aparicio, M. J, and Laroche S., (2011) An evaluation of the expression of the atmospheric refrac- 
tivity for GPS signals, JGR, 116, Dll, 16. 

Landau, L.D., E.M. Lifshitz (1987), The Classical Theory of Fields, Volume 2, Butterworth- 
Heinemann. 

Saastamoinen J, (1972) Contributions to the theory of atmospheric refraction, Bull Geod, 105:279- 
298 


6 







